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Abstract 

We reformulate the O(iV) sigma model as a loop model whose configura- 
tions are the all-order strong coupling graphs of the original model. The loop 
configurations are represented by a pointer list in the computer and a Monte 
Carlo update scheme is proposed. Sample simulations are reported and the 
method turns out to be similarly efficient as the reflection cluster method, 
but it has greater potential for systematic generalization to other lattice 
field theories. A variant action suggested by the method is also simulated 
and leads to a rather extreme demonstration of the concept of universality 
of the scaling or continuum limit. 
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1 Introduction 



The family of globally 0(N) invariant nonlinear sigma models, also called N- vector 
models, are very important statistical systems. For obvious reasons, in three space 
dimensions, they play a very prominent role in condensed matter physics. We here 
only mention the XY model (N = 2), relevant for the description of liquid helium, 
and the Heisenberg model (N = 3) for magnets. In high energy physics the four 
dimensional versions appear as effective field theories, for instance for pion physics, 
and a lot of interest focuses also on two (Euclidean) dimensions. This is motivated 
by the fact that these field theories are asymptotically free and share features with 
QCD like asymptotic freedom and dimensional transmutation with the nonpertur- 
bative generation of a scale like Aqcd- For example the study of a nonperturbative 
renormalized coupling constant in [TJ was an important preparation for the QCD 
Schrodinger functional methods [2J. 

With regard to the technique of Monte Carlo simulation, our main access to the 
models beyond perturbation theory, since about 20 years we are in an exceptional 
situation with regard to the O(N) models. The method of cluster updates [3] for 
embedded Z(2) (Ising) degrees of freedom [I] allows to painlessly enter the critical 
region, which unfortunately is in stark contrast to our possibilities in QCD. The 
hope that the cluster method would be widely generalizable was unfortunately 
disappointed in the following years, at least for high energy physics. In [5] even a 
kind of 'heuristic^] no-go theorem' was given concerning the generalization to sigma 
models with other spin manifolds. With the advent (or rather recognition^!) of [6] 
a completely different strategy to overcome slowing down has appeared: Simu- 
late the strong coupling graphs (to arbitrary order) instead of field configurations. 
Close to criticality the relevant field configurations are long-distance correlated. 
The cluster method, using an auxiliary percolation process, manages to execute 
collective moves that are thermodynamically appropriate for this case. Such moves 
seem to be difficult to find (and implement efficiently) in general!. The strong cou- 
pling graphs that are relevant at criticality are large and numerous. The clever 
idea of Prokof'ev and Svistunov [B] was to generate them not for the partition 
function alone, but to simultaneously consider the two point correlation. They 
have demonstrated in simple models that local deformations can pass between 
such graphs without significant obstructions and thus a relevant sample can be 
simulated. To avoid confusion we remind the reader of the following. In ordinary 
strong coupling expansions one takes the thermodynamic limit term by term and 
then the series usually has a finite radius of convergence which is often related 
to phase transitions. A finite lattice regularizes such singularities (for compact 

1 There are mathematical proofs which refer however to smooth field configurations. 
2 1 am indebted to Urs Wenger in this context. 
See [7] for a recent new proposal in this direction. 
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fields at least) and we can compute everything to in principle arbitrary precision 
by a convergent expansion for arbitrary couplings or temperature. While close to 
criticality this is impractical by conventional systematic expansion the stochastic 
evaluation is feasible. We may now also speak about an equivalence with another 
statistical system which incidentally has only discrete variables. Note that, al- 
though there are some similarities, this is not a complete duality transformation 
in the sense of Kramers and Wannier [8] . 

In a recently begun series of papers [9], [10] we have started to further work out 
the new approach. Beyond the Ising model we could apply it to fermions. Due to 
the sign problem this is at the moment still restricted to two dimensional systems 
like the Gross-Neveu model [TT]. A novelty that one has to appreciate is that the 
generated graphs can be adapted to the observables that one is interested in. One 
then needs several simulations for different quantities. We nonetheless see these 
dedicated simulations as a strength of the method. In [9], [12], [TT] it was found 
that this extra effort can result in enhanced precision for interesting observables. 

In this article we successfully extend the all-order strong coupling method to 
the class of O(N) nonlinear sigma models in arbitrary dimension. We first achieve 
this for the standard lattice action which allows to confirm our results by com- 
paring to other data in the literature. While the standard action has a relatively 
complicated all-order expansion we can define another action by insisting on a 
simpler expansion. Formally it can be argued to lie in the same universality class, 
but on the other hand it looks like a rather radical mutilation of the original spin 
model. We simulate its graphs and find that at least one universal result is repro- 
duced quite accurately and universality is confirmed. This flexibility will hopefully 
be useful to tackle further more complicated models in the future. 

In section 2 we develop the loop model equivalent to the all-order expansion of 
the O(iV) system. In section 3 we discuss how to represent the loop configurations 
in the computer and how to sample them. Extensive tests with the standard action 
are carried out in section 4. Section 5 discusses universality and the modified 
action followed by conclusions in section 6. In an appendix the limit N = 1 of our 
algorithm and its relation to previous Ising work is discussed. 

2 O(N) model as a loop ensemble 

We consider spin models with iV component spins s(x) of unit length located at 
the sites of a D dimensional hypertorus of length L M in the various directions. 
We refer all lengths to the isotropic lattice spacing thus putting a = 1. For the 
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standard lattice action the partition function with two field insertions reads 



Z(u,v) = J 



^i =(xy) s(x>s(y) s ( u y s ( v y 



The sum is over nearest neighbor links and the dots between pairs of spins mean 
O(N) invariant scalar products. The integrations employ the normalized O(N) 
invariant measure on the sphere, 



J dfi(s)f(s) = K N J d N s5(s 2 - l)f(s), K N ^ j dfi(s) = I. 



(2) 



Later on we shall need the corresponding single site generating function for a 
general source j a 



/oo 
d f ,(s)e^ = G N (j) = J24n;N}(j.j) 

n=0 



(3) 



which is essentially given by the modified Bessel function I N/2-1 and has expansion 
coefficients 

c\n- N] — r(iV/2) (4) 

The strong coupling expansion in (3 is generated by independently summing 
over an integer link field k(l) — 0, 1, 2, . . . , 00 in 



Z(u,v) = J2j 



Y[d»(s(z)) 



n (5) 



l=(xy) 



For a given configuration k the spin integral may now be written as 



X = 



d d 
dj a (u) dj a (v) 



n 

= {xy) 



d d 
djj(x) dj 7 (y) 



k(i) 



HG N (j(z))\^o. (6) 



We next introduce an auxiliary integer site field 



d(x) = 5 X , U + S X;V + ^2 k ( l ) 

l,dlBx 



(7) 



which counts the number of spins or respectively j-derivatives at x. It as well 
as X depends on u,v,k, of course, which we leave implicit for easier notation. 
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To produce a nonzero X, d(x) has to be even on all sites. Then the contribution 
becomes 



d d 



X> = ^TT^TT II TFTTTrTT 11^) -JW]^ 72 (8) 



where there are as many j factors as there are derivatives. In addition X' differs 
from X by dropping factors c[d(z)/2; N] for all sites z. The total number of terms 
in X' from taking all derivatives is 

M [u,v;k]=l[d(z)L (9) 

z 

The terms differ in their O(N) index contraction structures. To each of them 
there corresponds a graph A drawn on the lattice. There are k(l) lines between 
each nearest neighbor pair (xy) = I. At the 'interior of the sites' there is a kind 
of switch-board that sets up pairwise connections between all surrounding lines. 
Only at u and v two lines are left unpaired locally and instead are contracted 
with each other. Thus all lines are arranged in closed loops. The chain between 
u and v does not close geometrically (unless u = v) but closes with respect to 
O(N) contractions leading to a factor N as all other loops do. In the next section 
a visualization of such a graph or loop configuration will be given. Each graph 
represents a subset of the M.$ terms. This multiplicity is given by 

M[A] = -L (l[k(l)\J \[[d{x)/m dixV2 - (10) 

The last two factors correspond to permuting the pairs j-j at the sites and the two 
factors in each pair. In addition we consider the permutation of lines on the same 
link. For some graphs this leads however to an overcounting which is canceled by 
the symmetry factor S [A] . It is given by the number of elements of the group of line 
permutations which leave the connectivity of the graph unchanged. Apart from 
dealing with graphs embedded on a lattice it is similar to the symmetry factors 
that also appear for Feynman diagrams. We return to this issue in the discussion of 
our update scheme for graphs A. So far we have considered the loop configurations 
for given u,v,k. It is clear however, that the latter are also determined by the 
graph on the lattice. We may hence independently sum over graphs A G £ 2 which 
we define to include all possible locations u, v of the two 'defects' and all possible 
k(l) assignments to links that produce nonvanishing contributions. Then from (j5]) 
we generalize to 

z = J2p~\u-v)z(u,v)= ^Tp-V-^MA] (11) 

u,v Ae£2 
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where we have collected the whole loop weight into 



W[A) = N W M[A] 



n 

={xy) 



P S[A] 



]Jc[d(z)/2;N) 

z 

r(jy/2) 



-d(x)/2 



T(N/2 + d(x)/2) 



(12) 



Here u, v, k(l) and thus d(x) are now functions of A, and p is a positive weight to 
be chosen later. In the exponent |A| means the number of individual closed loops 
including the u-v chain. The loops in the configurations A that we sum over can 
overlap, intersect and backtrack. The weight depends on these features, the loops 
interact. In the next section we shall introduce an explicit parameterization of A 
together with an update scheme to simulate the loop model. 

As a by now standard next step [9] we introduce expectation values with respect 
to the new ensemble 



(13) 



Aec 2 



For the two point function of the original model there is the trivial relation 

5 aP Z(u,v) 



N Z(0) 



(14) 



where Z($) is the partition function without insertions (or u = v). This ratio can 
obviously be obtained from ( 1T3"j) as 



/ /n\ / \\ / \ \\^u—v,x)) 



where we have assumed the normalization 

p(0) = 1. 



(15) 



(16) 



It is convenient to in addition introduce expectation values referring to the subset 
of 'vacuum' configurations A G £2 that have u = v , 



«^(A)»o 



«<M) ' 

Such expectation values are independent of the choice of p. 



(17) 
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The internal energy density is equivalent to the average of the nearest neighbor 
correlation 



iV/ * — ' 

l={xy) 

where iVj is the total number of links. By differentiating Z($)=Z((5 UtV )) in both 
representations with respect to j3, the relation 



follows easily. Thus the average link occupation is bounded by (3. Numerical 
experience and large N considerations have shown that deep in the critical regime 
(3/N are typically numbers smaller than one. Thus, although the total length of 
the loops in A is unbounded in principle, as an extensive quantity it will in practice 
never exceed the number of links on the lattice (times N) by a large factor. We 
will come back to this point when we simulate. Another standard observable is 
the susceptibility 



which follows from contracting and summing over x in ()15p . 

The representation derived here allows to actually set A" also to non-integer 
values. The value N = 1 corresponds to the Ising model. The limit A" — > 1 of the 
present method does not immediately coincide with [B] , [S] . Therefore in appendix 
A we discuss the connection. 

3 Parameterizing and simulating the loop en- 
semble 

3.1 Parameterization: loops as lists 

It is often fruitful to first analyze graphs in an abstract manner and to separately 
consider their embedding on the lattice [13]. At the abstract level each graph 
A G £2 consists of lines and 2-vertices where two lines meet complemented with 
exactly two additional 1- vertices from the field insertions. We label all vertices with 
distinct integers. The lines are then naturally associated with pairs of integers. A 
valid embedding on the lattice associates lattice sites with vertices such that all 
lines map onto links, i.e. their index pairs refer to nearest neighbor sites. The 
1-vertices are at the sites u and v which can be anywhere on the lattice. 

We have found a representation for any embedded graph that we outline now. 
This is by no means unique. Our representation will actually be quite redundant 




(18) 




(19) 




(20) 



7 



by including extra information that will be found useful when we set up an update 
scheme in the next subsection. With any A we associate a list i which can be 
viewed as matrix Ey. It has one row for each vertex of the graph and there are five 
columns. In our convention the first and second row are permanently associated 
with the 1- vertices at u and v. The remaining rows then deal with 2- vertices and 
the row-indices % are taken as the graph theoretic labels of the vertices. In the first 
column in we encode the lattice site where the vertex is embedded. To this end 
we label the sites with integers 1,2, ... ,V in some order, with the lattice volume 
given by 

V = ]jL, (21) 
n 

Rows i of the list that are not in use for the given graph have in = 0. The 
configuration A in general holds many closed loops. The columns £i2 and £& are 
filled such they allow to travel around the loop containing a vertex % by following 
the pointers 

i -> i' = i a -> i" = i V2 f. (22) 

The loop passing through % can be traveled in two possible directions of which 
one is given by column 2. By using column 3 one obtains the other direction. 
Note that the loops of the O(N) models are physically unoriented, and we here 
encounter one of the redundancies mentioned before. The chain between u and v 
is treated analogously with the exception that the journeys are 1 — > i' = l\2 ~ * 

i" = U"i — > —^2 and 2 — > i' = £ 2 3 -» i" = — > -4- 1. From here on we 

call this special sequence of vertices and lines the active loop@ with the remaining 
ones being called passive. Column 4 just holds a flag whose values distinguish 
vertices in passive loops (one) from those in the active loop (zero). Finally column 
5 is arranged such that the following problem can be solved efficiently, i.e. without 
searching the whole list: for a given lattice site^l x find all vertices (row-indices) 
embedded at this site for the present A. With an additional entry-list e(x) the 
problem is solved by following the sequence 

x _> e {x) = i ->■ i' = i ib -> i" = i Vh (23) 

This chain ends when £j$ = is encountered. During the later update a vertex 
can be removed from a graph. If it corresponds to row i, we set in = in such a 
case. Such a row can continue however to function in (I23p . As a consequence, in 
0231) we can step through such lines which are still needed even without holding 
a vertex. In addition to I and e we also store the auxiliary field d(x), although it 
could be constructed from £. 

4 Remember it is a closed loop in the O(N) sense yielding a factor N. 

5 By x we designate both the geometric location as well as the counting label given to it. 
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Figure 1: A small loop configuration to illustrate its representation by a list. 



We now discuss as an example a 3 x 3 lattice taken from an 0(3) run. It is tiny 
to save space but it actually demonstrates most features. A loop configuration is 
shown in figure [TJ We have displaced the loops in the third direction to disentangle 
overlapping loops. On the 'roof we show the labeling of lattice sites, the other 
integers are vertex labels. The solid line, connecting 1 with 2 (via 7 and 11) is the 
active loop and there are two additional closed passive loops which overlap at site 
5. 

As already mentioned there is no hard bound on the number of rows/vertices 
that are needed in t. In addition, as A is updated (modified), vertices are added 
and eliminated and thus rows of I are freed (possibly with the exception of the 
entry in column 5) and new ones are required. As the available storage is finite, 
this requires some management which may seem difficult at first sight. Luckily 
this problem can be handled rather easily. 

As argued before the average total number of vertices will be of order NxNi. As 
an extensive quantity its fluctuations are found to be only of the order y/N x iVj. 
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We thus found that reserving space for 



{£ i , j=l ,..., 5 },i = l,2,...J £ xNN l (24) 

lines with fi of order unity leads to completely negligible probabilities to ever 
exhaust this space in any feasible simulation. These storage requirements are 
quite similar to ordinary simulations. By observing the fluctuations of J2 X d(x) 
one can easily demonstrate failure probabilities like for example 10~ 1000 which can 
be tolerated. We come back to this in section 4. 

For recycling list entries we keep another list with a reservoir of indices of 
completely unused rows that are available for new vertices. There is a subtlety 
here. As mentioned before a line not carrying a vertex anymore can still be relevant 
as a 'stepping stone' in (1231) . Such rows we call 'unused' as opposed to 'completely 
unused'. If we run out of completely unused rows during a simulation we start a 
recycling routine which runs through the whole list. In this process all £^ where 
li\ ^ and the associated entries in e are recomputed and all rows with £n = 
acquire the 'completely unused' status. It turns out that the time spent for list 
recycling is a negligible fraction of the total in practice. We end this subsection by 
reproducing in table [T] the lists £ and e associated with the configuration in figure 

m 

1 1 7 

2 4 

3 21 

4 5 24 

5 6 10 

6 9 

7 4 11 

8 22 

9 8 15 

123456789 
12 3 7 4 5 11 9 10 
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1 





24 


24 


1 


15 


15 


1 


8 


1 








1 





16 


1 








10 


1 






10 


9 


9 


5 


1 


22 


11 


7 


2 


7 





17 


12 





2 


1 








15 


5 


5 


9 


1 


6 


16 





17 


17 








17 





7 


7 








22 





9 


5 








24 


2 


4 


4 


1 






Table 1: The list £ corresponding to the configuration in figured] (upper two parts). 
It has been augmented by the leftmost (zeroth) column exhibiting the row indices. 
Completely unused lines have been omitted. The lower list is the entry table e(x). 
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3.2 Monte Carlo algorithm 

We now propose an algorithm to simulate the 0(N) loop ensemble (IT3|) . We define 
a number of separate update steps such that each of them fulfills detailed balance. 
They will then be iterated in some order as the final update procedure. The 
moves are all Metropolis proposals for which we quote the ratio q which controls 
the acceptance probability min(l, q) in each case. We need to introduce the notion 
that the active loop between u and v is called trivial if it contains no 2-vertex and 
also u = v coincide. 

I. Extension and retraction: We choose with equal probability between 2D + 
1 possible proposals to move u by one lattice spacing with a concurrent 
adjustment of the active loop. In the first 2D cases it is extended with u 
moving to one of its neighbors u with the amplitude ratio 

qcxt ~ N + d[u]p(u-vY [ ' 

In the last case u is retracted by one link along the active loop with the ratio 

N + d[u]-2p(u-v) 

No move is made in the last case if the active loop is trivial. 

II. Re-route: We here want to change the O(N) contraction or line connectivity 
structure at u. We have to distinguish a few cases. In the cases not covered 
below no move is made. 

i. The active loop is trivial and d(u) > 2. In this case the we pick a 2-vertex 
at u and replace it by the two 1-vertices, assigning u and v randomly to 
the two lines. The acceptance ratio is 

*. = ^ (27) 

Note that upon acceptance the active loop becomes non-trivial and the 
loop number |A| is reduced by one. 

ii. The active loop is not trivial but u = v holds. In this case we make the 
move inverse to i with ratio (g r er) _1 - 

iii. We have u ^ v and d{u) > 2. We pick with equal probability one of the 
lines connected to any of the 2-vertices at u and propose to redirect it 
to the 1-vertex. The line previously connected to the latter is rewired to 
the newly created 'opening' at the 2-vertex. For the acceptance decision 
we need to distinguish further sub-cases. In figure [2] we give a hopefully 
helpful illustration of the various moves. 
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a) The chosen 2-vertex belongs to a passive loop. The latter then gets 
inserted into the active loop, |A| is reduced by one, and the ratio is 
l/N in this case. 

b) The chosen 2-vertex belongs to the active loop, which self-intersects 
at u, and the chosen line leads towards the 1-vertex at v. In this case 
a new passive loop is detached and |A| goes up by one. The ratio is 
N. 

c) The chosen 2-vertex belongs to the active loop and the chosen line 
does not lead to the 1-vertex at v . In this case the active loop is just 
re-ordered and the ratio is one. 

III. Kick: Do nothing unless the active loop is trivial. If this is the case we pick a 
random site x and propose to move both u — v to x. The active loop remains 
trivial. The acceptance amplitude ratio is 

_ N + d[u]-2 
qklck ~ N + d[x] ■ (28) 




Figure 2: Illustration for the re-routing moves. The dashed line encircles elements 
associated with one site. The x stand for the 1-vertices at u, v while a line both 
entering and leaving the dashed circle represents a 2-vertex. There could be more 
'spectator' 2-vertices at the site which are not drawn for clarity. 
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Figure 3: A configuration which locally contributes a factor 3! x 8 to the symmetry 
factor S. Vertices at u are those inside the dashed circle. Under moves of type 
Iliiia S can be reduced by factors 1, 3 or 4 in this case. 



The symmetry factor S [A] in (|T2"|) can only change during steps II, because the 
active loop is distinguished and its changes do not influence the symmetry of the 
graph A. If in II the symmetry changes, this is effectively taken into account 
by asymmetric proposal probabilities for the forward and backward process. We 
explain this by an example given in figure |3l In the figure contributions to the sym- 
metry factor, where line permutations leave the graph unchanged, are indicated. 
In step Iliii we first pick a line from one of the 2-vertices to then connect it to 
the 1-vertex at u. If lines participate in the symmetries several choices lead to the 
same proposal A — > A'. More precisely the new graph is reached from <S[A]/«S[A'] 
such choices. In this way the total transition probability in Iliiia is 



In this way we have detailed balance with respect to (I12p including S. Analogous 
considerations apply to the other moves in step II. 

It should be easy to imagine now how the available information in the list i 
described before is useful during the update steps. For instance the active/passive 
flag helps to discriminate between the sub-cases of Iliii based on locally available 
information. It is also clear that precisely analogous steps can be defined around 
v instead of u chosen above. After some brief experiments we have arranged the 
updates in the following way 




(29) 



while the reverse process Iliiib proceeds with 




(30) 



1 Iteration := (LjILjLjIIy III) 



NxV/2 



(31) 
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This was partly based on aesthetics and symmetry. For instance making only u 
steps or dropping III made little difference for our observable^). All Metropolis 
acceptance rates are well above 0.5 with the exception of the extensions step in I. 
The latter is close to (2D)~ 1 . For simulations with the modified action in section 
5 it will rise to 1 however. 

The mathematical proof of ergodicity is the usual one: with a nonzero proba- 
bility any configuration A can be transformed to the trivial empty one and then 
evolved to any other A'. It is also the empty lattice from which we start all simula- 
tions. Another observation is the following. A correct algorithm is also given with 
only the steps i and ii contained in II, thus re-routing only for u = v. We have 
indeed confirmed some correct results with such an update, but it is accompanied 
by severe critical slowing down. We shall see that this is completely eliminated by 
the additional stepsl- 




Observables are accumulated after each l u and I„ step during the update. In 
most runs these contributions are stored separately for each iteration, and later 
an off-line autocorrelation analysis [2] is carried out for these time series. We 
then arrive at integrated autocorrelation times in units of iterations fl3TT) . During 
the mostly local steps needed to update the list £ we sometimes have to trace a 
closed loop, for example to adjust the active/passive flag. Because of this it is not 
obvious at this point that the effort for one iteration scales strictly proportional 
to the lattice size V. We shall come back to this point in the next section. 

4 Numerical experiments with the standard ac- 
tion 

We have first made series of runs for the loop model representing the standard 
lattice action discussed up to here. Most of these results can hence be directly 
compared to numbers in the literature. In table [2] we have listed run parameters 
and the observed integrated autocorrelation times for a few observables. For the A- 
series the (3 values of [T5] have been adopted together with lattice sizes to arrange 
for mL « 8 to hold. According to [T6] our massgap then differs from the infinite 
volume one only at the level of which will be below our errors. The weight 
p in (fl3|) is chosen to roughly anticipate the decay of the two point function. 
The relative error of the true correlation (fl"5l) is then constant or even shrinking 
as we explore the exponential fall-off over a long range. We refer to the detailed 

6 If both is done, v does not move anymore. Observables where we average over translations 
still seemed to assume correct values. 

In the zoological interpretation as a worm algorithm [5] step I corresponds to the normal 
development of a worm. II deals with the asexual reproduction by detaching parts of its body 
(Iliiib). Sadly, the O(N) worm sometimes also devours its offspring with probability 1/N (Iliiia). 
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discussion in [9] which applies here unaltered. In this paper we use the free massive 
lattice propagator to determine p, 

(-A + M 2 )}\x) = 4,o p{x) = f{x)/f{0), (32) 



where A is the standard nearest neighbor lattice Laplacian. Using the fast Fourier 
transform on one lattice direction after another, its construction costs negligible 
0[DL D \nD] operations. Thus the column for M completely determines p(x). 



run 


N 


D 


L 





M 


T~mt,K 


T int,x 


7~irit,m 


Tint,|A| 


CPU 


Ai 


3 


2 


56 


1.4 


8L- 1 


1.84(3) 


0.518(3) 


0.678(6) 


1.46(2) 


1.50 


A 2 


3 


2 


88 


1.5 


8L- 1 


1.74(3) 


0.515(3) 


0.691(6) 


1.39(2) 


1.58 


A 3 


3 


2 


152 


1.6 


8L- 1 


1.72(3) 


0.515(3) 


0.698(6) 


1.38(2) 


1.75 


A 4 


3 


2 


276 


1.7 


8L- 1 


1.71(3) 


0.514(3) 


0.702(6) 


1.37(2) 


2.01 


A 5 


3 


2 


518 


1.8 


8L- 1 


1.65(3) 


0.510(3) 


0.689(6) 


1.34(2) 


2.56 


B 


8 


2 


128 


5.2 


8L- 1 


1.94(3) 


0.506(3) 


0.571(4) 


1.48(2) 


2.51 


Ci 


1 


2 


128 


/r 


p=l 


2.42(6) 


0.927(14) 




2.41(6) 


1.42 


c 2 


1 


2 


256 




p=l 


2.80(7) 


1.005(16) 




2.83(7) 


1.49 


Di 


3 


3 


32 


pF 


p=l 


4.08(9) 


0.818(8) 




3.02(6) 


1.50 


D 2 


3 


3 


64 




p=l 


4.92(12) 


0.910(10) 




3.63(8) 


1.83 


Ei 


3 


2 


16 


1.779 


2L- 1 


0.824(8) 


0.536(4) 


0.776(8) 


0.72(1) 


1.95 


E 2 


3 


2 


32 


1.779 


2L- 1 


0.900(10) 


0.527(3) 


0.842(9) 


0.79(1) 


2.21 



Table 2: Run parameters and autocorrelation times. The statistics for each run 
consists of 10 6 iterations. For A...D the geometry is symmetric L M = L, while in 
E it is elongated in 'time', L = QLi = QL. The critical value pc — = 0.693002 was 
taken from [T7] . 



In table [3] the corresponding mean values and errors are compiled. They refer 
to the mean bond occupation (fl9l) . the susceptibility (1201 . and the loop number 
appearing in (fl~2l) . In addition we have recorded the two point function 

G{t) = ((p(u - v)[5 t ,u -v + St^-vA))- (33) 

In this formula, valid for the D = 2 symmetric geometry, we sum over both 
directions and 5 is taken L-periodic. In addition we average over reflections. We 
have checked that these contributions are not completely, but largely statistically 
independent. 

In figure H] the errorbars show effective masses derived from the ratio of suc- 
cessive time-slice correlations by matching with cosh(m e fr(t — L/2)). We see the 
expected long plateau and it seems completely safe from excited states errors to 



15 



run 


K 


X 


™— i 
m 


//I A l\\ T/— 1 

((|A|))o x 1/ 


t > r 

Her. 


Ai 


0.78701(8) 


78.75(16) 


6.876(5) 


0.34590(4) 


[15j 


A 2 


0.90246(6) 


175.38(41) 


11.053(8) 


0.34818(3) 


[15J 


A 3 


1.01715(4) 


448.3(1.1) 


19.030(13) 


0.34604(2) 


[15j 


A 4 


1.12920(2) 


1269.4(3.5) 


34.530(25) 


0.34244(1) 


[15J 


A 5 


1.23829(1) 


3855(11) 


64.872(66) 


0.33927(1) 


[15j 


B 


3.38678(6) 


431.54(80) 


18.096(9) 


0.88431(3) 


1 -i r-i 1 

[18J 


Ci 


0.31264(8) 


1.095(7)L 7 / 4 




0.13190(4) 




c 2 


0.31217(5) 


1.104(8)L 7 / 4 




0.13216(3) 


i 


D x 


0.23015(2) 


1.1055(19)L 2 




0.19294(1) 




D 2 


0.22910(1) 


1.0771(17)L 2 




0.193582(4) 




Ei 


1.2178(2) 


242.1(1.4) 


15.115(26) 


0.33892(11) 


[19] 


E 2 


1.21662(9) 


632.5(3.1) 


25.115(37) 


0.33867(5) 


[19J 



Table 3: Values of some physical observables defined in the text. In the last column 
we cite references, where data consistent with the ones here can be found. 



finally extract the mass from a fit, for which the horizontal line shows the range 
and the value in the form of an drier error band. For the fit we have minimized 
over the shown range the function 

^ [G(t)-ccosh(m(t-L/2))] 2 

5GW ( ) 

with respect to c and m. More precisely, we first determine the error SG(t) by 
analyzing G(t). As expected SG(t)/G(t) hardly grows with t. Then these errors 
are used to define via the minimization of ( I34j) m as a function of the primary 
correlation data. An error for this derived observable is estimated as discussed 
in [H]. These are the values quoted for m in table |3j Although the effective 
masses fluctuate around the plateau (see figure HJ), the G{t) cannot be expected 
to be completely uncorrelated at neighboring t values. While fl3"4"|) is the so-called 
uncorrelated X 2 we actually see values between 0.1 and 0.4 per degree of freedom. 
This is not extremely small, we collect much more independent information than is 
customary in standard simulations. In fact the error of the fitted mass is between 
2 (for Ai) and 4 (for A 5 ) times smaller than that of the effective mass at the 
beginning of the fit range alone. We have found that the mass values and their 
errors change only within their errors, if we replace the weight 5G~ 2 in (154]) by one 
flat in t. Plots similar to figure H] arise for all lattices where a mass is quoted. The 
data in [15] allow to roughly estimate which statistics was invested for the errors 
quoted in units of steps per spin. We conclude that for the estimation of the mass 
gap, the present method is quite competitive with the reflection cluster algorithm 
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0.5 1 1.5 2 2.5 3 3.5 4 

m t 

Figure 4: Effective masses and fitted mass (band between the lines) from run A4 
(L = 276). To avoid cluttering only every second m e ff is shown. 



[I] with improved estimator [15]. This is not quite so for Xi which could possibly 
profit from a different choice of p. 

In the series C we investigate the exactly solved two dimensional Ising model 
by simply setting N — 1 in our loop code. The efficiency is in fact quite similar 
to the simulations in [9] although the sampling of the contraction structures in 
({TBI is an in principle unnecessary complication as iV' A ' equals unity in this case, 
see the appendix for further remarks. In the series D we take the 0(3) model to 
three dimensions and simulate at the critical point = 0.693002 determined in 
[17] . Beside the two runs quoted we have also reproduced data from [20]. The C 
and D cases demonstrate the very mild critical slowing down of our simulations at 
criticality where L is the only scale. 

With E we turn to studies of the renormalized coupling 

g 2 = m(L)L (N = 3) (35) 

that has been introduced in pQ. Here m(L) is the mass gap of the transfer matrix 
of spatial size L. In [TJ it was argued that at least in the perturbative regime free 
boundary conditions in the time direction help to isolate the first excited state 
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Figure 5: Similar to figure IU but for simulation E 2 (L = 192, Li = 32). 



from the rest of the spectrum. The loop model could easily be modified to these 
boundary conditions. We here prefer however to make the time direction long 
(Lq = QL) to achieve the isolation. This is facilitated again by our small errors at 
large separation. An advantage of this approach is that exact translation invariance 
is kept for both directions. In figure [5J, analogous to HI we demonstrate how the 
gap is extracted also in this case. Here we took M > m and correspondingly 
'oversample' [9] large distances with errors shrinking with growing t (except very 
close to Lq/2). Recently in [21] and [19] a thorough study of cutoff effects in 
the step scaling function for (1551) was made. In this context very accurate data 
were produced including our lattices E 1; E 2 where we obtain g 2 = 1.0586(18) and 
~g 2 = 1.2741(19). In this case a special estimator (22] is available in the spin 
formulation such that our accuracy achieved here cannot really compete with [19], 
but the more accurate values are consistent with ours. 

We now make some general observations on the loop configurations observed 
in our runs. 

An exceedingly encouraging observation on our compiled autocorrelation times 
is, that in units of iterations there seems to be almost no critical slowing down 
for any of the series and quantities studied here. The typical bond occupations 
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K resemble E, see ( 1T9|) . and are functions of /3 with only a weak dependence 
on Lfj_. They are larger in lower dimension. The number of loops |A| is strictly 
proportional to V and also similar to E. Both quantities grow with N for similar 
correlation lengths as one would expect. 

We finally return to the question of the computational complexity of one iter- 
ation. Here the last column 'CPU' in tabled is of interest. It gives the execution 
time in \i sec of one micro step, i.e. the time for one iteration divided by NV. This 
is of course a highly non-universal implementation and processor dependent quan- 
tity of not much interest. It is quoted however, because its relative change between 
the runs may be of some more interest. A constant value for CPU would suggest 
a scaling behavior like for sweeps of local algorithms. The additional growth that 
we see within the simulation series represents a small effective slowing down in 
CPU time units. In the A series we really take a scaling limit. From Ai to A5 the 
correlation length changes by a factor 9.4. The extra CPU time factor is 1.7. The 
reason for this extra growth indeed comes from the 'looping' steps as discussed 
at the end of the last section. For a discussion of the true asymptotic dynami- 
cal behavior, where this component will probably eventually dominate, one could 
investigate the distribution and scaling behavior of the perimeters of the loops 
in analogy to percolation cluster sizes. We however do not try to determine dy- 
namical exponents. We just conclude that the new type of simulations is efficient 
enough to generate all data cited here in a few hundred hours on a PC with a code 
that still allows for ample speedup. It is not clear at the moment if a modified 
list parameterization and/or algorithm could avoid this more than linear with the 
lattice volume growth of the CPU time per iteration. Note that we have been dis- 
cussing the continuum limit throughout. For the thermodynamic limit (V — > 00 
at fixed (3) strict linearity is expected for the massive theory. 

We close by a few remarks about our specific implementation that is reflected 
in the quantity CPU. Our code is written in and running under mat lab. Only 
the random numbers are imported from the C-code [23] and we use luxury level 
two throughout. We used up-to-date (2008) PCs with four cores that we employ 
for trivial parallelization thus speeding up the time following from table [2] by a 
factor four. We hence always have four replica which allow for another reassuring 
consistency check on the error determination by always monitoring the Q-values 
[H]. It would be clearly possible to speed up our runs by a large factor by writing a 
dedicated C-code. In particular the data type 'pointer' seems ideal for the handling 
of loops by lists. Like for cluster simulations, the present algorithm presumably 
does not lend itself very naturally to a nontrivial parallelization. 
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5 Simplified actions 



For the 0(N) model we have succeeded to simulate the untruncated strong cou- 
pling expansion for the standard nearest neighbor lattice action. This could be 
achieved because the weight needed to integrate out the spins ([3]) is relatively sim- 
ple in this case. We anticipate that in other models that we shall want to treat 
similarly this may be more involved. It would be easier to handle the graphs and 
to tabulate the required weights if we could limit the overlapping of loops. For the 
O(N) model this would amount to constraining the sum in (j5J) and in the formulae 
following from it by 

k(l) = 0,1,2,... ,k max . (36) 

In the simulations reported above the single link occupations k(l) exhibit a Poisson 
type distribution. Hence, for the mean values given, they never get much beyond 
about 10 and a /c max of this size would have very little effect. From here on we 
shall however investigate the most radical possibility by setting fc max = 1. Then 
we effectively replace (P) by 

Z(u,v)= ! Y[dfi(s(z)) s(u)-s(v) l[ [l + Ps(x)-s(y)j (37) 

J I Z J l = ( X y) 

and use a tilde for quantities which refer to this action. It is clearly still ultralocal 
and, if such a system becomes critical at all, we would expect to be in the same 
universality class as before. The Boltzmann factor is strictly positive for \(3\ < 1 
only. In the Ising model at N = 1 the truncation fc max = 1 is related to expanding 
in (3 = tanh (3 instead of (3. Then, with this identification, we compute exactly the 
same correlations on every finite lattice. 

A weight like ( 1371) is actually not new in the literature. In [21] this action was 
put on a honeycomb lattice in two dimensions. Since only three links meet at a 
site there, the strong coupling graphs simplify as they cannot overlap or intersect, 
similarly to the Majorana fermions in [10]. In this way relations with certain 
discrete models can be derived. This is elaborated in [25] where critical indices 
are derived for —2 ^ N ^ 2. In both publications it is conjectured that, although 
unusual, the Boltzmann weight in ( 13T|) should lead to the O(iV) universality class. 
We here check this for the 0(3) model in D = 2 (on our usual square lattice) by 
computing the step scaling function (SSF) of pp. 

It involves pairs of lattices 

S(2, u, L- 1 ) = m(2L)2L\ m(L)L=u (38) 

where the side condition determines the value (3 (or (3) to be used on both lattices. 
For each u the SSF is expected to have a universal continuum limit 

<r(2, u) = E(2, it, 0) = lim E(2, u, L' 1 ) (39) 

L— >oo 
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where a is equivalent to the Callan Symanzik beta function for the renormal- 
ization scheme defined by the coupling (155)1 . It is quite amazing that by Bethe 
Ansatz techniques the continuum SSF a could be computed exactly [26J. For our 
universality check we pick the particularly popular point p], [19] 

wo = 1.0595, ct(2, u ) = 1.261210. (40) 

The exact result has been confirmed with extreme numerical precision in [TP"] . 

A few short experiments immediately showed that > 1 is required to achieve 
any sizable correlation length in lattice units. Then the weight in ( 157)1 oscillates 
and presumably has a severe 'sign problem'. The loop partition function Z which 
is just as pip but omitting terms with any k(l) > 1 continues to be a sum of 
positive terms, of course. It is also plausible that loops in A grow with rising 
and thus encode correlations of growing range. The counterpart of (1T91 now reads 

In the scaling region we shall find mean bond occupations not far from 1/2. Due 
to the non-positive weight it is however not possible to draw any immediate con- 
clusions on the nature of 'typical' field configurations in the spin representation. It 
also seems hard to imagine to set up a bare perturbative expansion on the lattice 
with this action. This is perhaps somewhat reminiscent of the constraint model 
in [27] where the Boltzmann factor acts ferromagnetically only in the form of a 
step function on the angle between neighboring spins. If features of the long range 
physics are nevertheless described by renormalized perturbation theory, this must 
be viewed more like an effective field theory matched to these lattice models. 

We now come to our concrete datc@ taken allowing k(l) = 0,1 only. For the SSF 
we have to first determine the values that yield ~g 2 =Ho on a series of lattices. Such 
results are listed in table HI We have measured the derivative dg 2 /d(3 and used 
it so solve the problem of tuning (6th column). Our simulations are so close to 
the target that further terms in a Taylor expansion and the errors of the measured 
derivative play no role. We note that is larger than the corresponding in E 12 
and that there are about four times fewer loops now. The effective mass plots were 
again inspected. They all look qualitatively indistinguishable from figure Again 
M = 2/L determined p{u — v) for us. 

In table [5] we list the corresponding doubled lattices required for the SSF. The 
error quoted for E(2, u , 1/L) combines the statistical errors of both the small and 
the doubled lattice. 

8 Due to a programming error the data in this section had to be revised in comparison to the 
previous arXiv versions. 
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J 

L 


p 


K 


/I a l\ \r— l 
(|A|)o x 1/ 


tt2 
9 


P{Uq) 


IrU 


6 


1.7851 


0.41112(11) 


0.09098(5) 


1.0622(9) 


1.7904(18) 


1.42 


8 


1.9107 


0.43875(9) 


0.08468(5) 


1.0612(10) 


1.9144(21) 


1.46 


12 


2.1114 


0.47897(6) 


0.08125(4) 


1.0622(10) 


2.1185(26) 


1.54 


16 


2.2841 


0.51044(5) 


0.08060(3) 


1.0636(10) 


2.2955(28) 


1.74 


24 


2.5742 


0.55930(3) 


0.08073(2) 


1.0590(10) 


2.5724(36) 


1.93 


50 


3.1055 


0.64402(2) 


0.08165(1) 


1.0597(11) 


3.1063(32) 


2.72 



Table 4: Lattice data with k max= i for the 0(3) model in D = 2 with Lq = 6L\ = 
6L. Each line represents 3 x 10 6 iterations. 
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(3 
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E(2,«o,l/L) 


CPU 


12 


1.7904 


0.40987(10) 


0.07897(5) 


1.2300(17) 


1.2300(22) 


1.48 


16 


1.9144 


0.43832(7) 


0.07880(4) 


1.2304(17) 


1.2304(23) 


1.53 


24 


2.1185 


0.47978(5) 


0.07913(3) 


1.2424(17) 


1.2424(23) 


1.63 


32 


2.2955 


0.51219(4) 


0.07966(2) 


1.2448(18) 


1.2448(23) 


1.93 


48 


2.5724 


0.55882(3) 


0.08045(2) 


1.2547(18) 


1.2547(25) 


2.15 


100 


3.1055 


0.64396(2) 


0.08162(1) 


1.2606(19) 


1.2600(30) 


3.34 



Table 5: The doubled lattices of table HJ 10 6 iterations each. 



In figure [6] we plot our SSF data. The left panel shows £ versus L~ 2 with 
the star at zero giving the exact continuum limit It is very obvious that 

the data deriving from the action in ( 13 7p 'know' about the universal value. In 
[19] a detailed investigation on the Symanzik effective field theory for the lattice 
artifacts in the two dimensional O(N) models was carried out. The result is that 
the leading term for 0(3) is expected to be of the form In 3 L/L 2 with sub leading 
contributions having smaller powers of the logarithm. This has motivated us to 
plot (S — cr)L 2 / In 3 L versus In" 1 L. The two lines are linear fits with acceptable 
X 2 values. The dashed line differs from the solid one by including/omitting the 
coarsest lattice pair. We see how plotting matters: what is a short extrapolation 
in the left plot is an uncomfortably long one in the right panel. But the artifacts 
for our action are clearly compatible with the theoretical form. In addition it is 
reassuring to note that any non-absurd extrapolation in the left plot would not miss 
the continuum value by much. The finest lattice pair agrees with the continuum 
value within our small statistical error. 
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Figure 6: Data for the step scaling function with the k max = 1 action. The star in 
the left plot is the known exact answer. The right plot is relevant for the theory 
of scaling violations a la Symanzik in 



6 Conclusions and outlook 

We have successfully extended the method of simulating strong coupling graphs 
instead of fields to the O(N) models. They are in this way first reformulated as 
a loop model which is then simulated. The resulting setup of reformulation plus 
algorithm is, at least for the mass gap as an observable, similarly efficient as the 
reflection cluster method. A variant of the lattice Boltzmann factor, suggested 
by simplifying the loop ensemble, was simulated and demonstrated to yield the 
same continuum step scaling function for a particular value of the coupling. Hence 
universality is confirmed here in an interesting case. A direct simulation as a spin 
model would presumably be difficult due to sign oscillations of the Boltzmann 
weight. 

The main goal of this work was not primarily to have another simulation 
method for the O(N) models. The main motivation was that the further ex- 
tension of this technique is possible in a more systematic way than for the cluster 
method. After all, many systems possess a strong coupling expansion. The O(N) 
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model is simply a very well studied case which has taught us important lessons 
about the general approach. As a byproduct one could now however study the 
O(N) models for non-integer N and perhaps one could also consider taking the 
N — > oo limit (at fixed (3/N) of the algorithm similarly to the time continuum 
limit in [28J. Also, as already indicated above, the study of the individual loop 
distribution in the spirit of percolation theory could be of some interest. 

A modest next step in the main line of the project will be to consider CP(iV— 1) 
spin models, for which no efficient cluster algorithm exists, at least for the classical 
models (see however [29]). A more long distance goal remains the simulation of 
lattice gauge theory as an equivalent surface model. 
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A The Ising limit N = 1 



For N — 1 the O(N) system becomes the Ising model. In this case the function 
([3]) is Gi(j) = cosh(j) and hence 



c(n: 1) 



1 



such that 



(2n)! 

]Jc[d(z)/2;l} = ^ 



Mo 



may be used in (FI2"]) . With no more dependence on |A| 
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Aec 2 



M [u,v; k] 



n 



may be rewritten as a sum over u, v, k 



* = £/>-'< 



u — v) 



u,v,k 
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ip 



d(x),even- 



(42) 
(43) 



(44) 



(45) 



which is the expected /3 expansion for the Ising model. Note that u, v enter via 
the relation (JTj). While the form fTPTl) contains the sum over different contractions 
it is also correct with N = 1. All contractions are equally likely in this case, the 
sum over them factorizes off. 
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It is interesting to interpret our general update from the point of view of the 
ensemble ( l4"5j) . To obtain the effective transition probability from step I we average 
over possible contractions at u and sum over the resulting ones at u. This leads 
(for p = 1) to the following update sequence 

• With probability 1 - l/(2£> + 1): 

— Pick a link I = (uu) around u in one of the 2D directions with equal 
probability. 

— Move u — y u [and k{l) — > k(l) + 1)] with probability 
min[l,/3/(d(«) + l)]. 

• With probability 1/(2.0 + 1): 

— Pick a link I = (uu) around u in one of the 2D directions with probabil- 
ity k{l)/(d{u) — 1). These probabilities add up to unity for the 2D link 
directions if u ^ v holds. For u = v there is a probability l/(d(u) — 1) 
that no direction is picked. In this case no move is made. 

— Otherwise move u — » u [and k(l) — y k(l) — 1] with probability 
min[l, (d(u) - l)/0] 

The ratio in the third sub-item is the fraction of all possible pairings that connects 
u to a link in the given direction. It is not hard to see that these implied transitions 
for u, v, k obey detailed balance with respect to (14"5"1) . 
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